{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению. Сессия № 3\n",
    "\n",
    "### <center> Автор материала: Латышев Евгений (@EvgeniyLatyshev)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center> Индивидуальный проект по анализу данных </center>\n",
    "### <center> Прогнозирование уровня производства конфет в США </center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = 20, 15\n",
    "import seaborn as sns\n",
    "from scipy import stats\n",
    "import pylab\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import statsmodels.api as sm\n",
    "from itertools import product\n",
    "import fbprophet\n",
    "from tsfresh import extract_features, select_features, extract_relevant_features\n",
    "from tsfresh.utilities.dataframe_functions import impute, make_forecasting_frame\n",
    "import xgboost as xgb\n",
    "from sklearn.metrics import roc_auc_score, mean_absolute_error\n",
    "from hyperopt import hp, fmin, tpe, STATUS_OK, Trials\n",
    "\n",
    "RANDOM_STATE=42"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Описание набора данных и признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Каждый год люди съедают огромное количество конфет, которое сначала необходимо произвести. Очевидно, что уровень производства постоянно меняется. Например, в канун нового года или на хеллоуин съедается гораздо больше сладостей, чем в другие времена года, что вносит элемент сезонности. Также, в долгосрочной перспективе, уровень производства непреклонно растёт вместе с ростом населения, что говорит об имеющейся трендовой составляющей. Эти свойства делают рассматриваемый показатель отличным подопытным с точки зрения анализа временных рядов.\n",
    "\n",
    "Датасет, полученный с [Kaggle](https://www.kaggle.com/rtatman/us-candy-production-by-month) содержит данные об уровне производства конфет в США с января 1972 по август 2017 года.\n",
    "\n",
    "Каждая запись характеризуется двумя значениями:\n",
    "- observation_date - дата в формате YYYY-MM-DD\n",
    "- IPG3113N - [индустриальный продуктовый индекс](https://fred.stlouisfed.org/series/INDPRO). Это универсальный индекс уровня производства, который измеряется как % от уровня производста 2012 года.\n",
    "\n",
    "В данной работе рассматривается задача прогнозирования будущего уровня производства конфет по имеющимся предыдущим наблюдениям. В качестве горизонта прогнозирования выбран временной интервал в 24 месяца."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Первичный анализ признаков, визуальный анализ признаков, поиск зависимостей"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('../../data/candy_production.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сперва приведём данные к более удобному виду."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.index = pd.to_datetime(df['observation_date'], format='%Y-%m-%d')\n",
    "df['value'] = df['IPG3113N']\n",
    "df = df.drop(['observation_date', 'IPG3113N'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь проверим данные на наличие пропусков и выбросов."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.isnull().sum(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видим, датасет содержит 548 значений в промежутке от 50 до 139. Судя по квантилям, данные относительно чистые и не требуют предобработки."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Далее можно выполнить проверку на нормальность и скошенность. Используем критерий Шапиро-Уилка, Q-Q plot и проверку на скошенность из библиотеки scipy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.shapiro(df['value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.probplot(df['value'], dist=\"norm\", plot=pylab)\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['value'].plot(kind='hist', bins=25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.skewtest(df['value'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Распределение целевой переменной близко к нормальному распределению, но немного скошено."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Пора посмотреть на график самого временного ряда."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.tsplot(data=df['value'], time=df.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим, что нам даст группировка наблюдений по месяцам и годам."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['year'] = df.index.year\n",
    "df['month'] = df.index.month"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.boxplot(x='year', y='value', data=df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.boxplot(x='month', y='value', data=df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Не считая спада 1975 года, количесво произведённых конфет продолжало уверенно расти вплоть до начала нового тысечелетия, что, возможно связано с ростом населения США, которое за эти годы увеличилось с 210 млн. до 282 млн. Однако, начиная с 2000 года наблюдается некая стагнация, выраженная сравнительно небольшим разбросом производства в эти года, а с 2005 года уровень производства и вовсе начинает падать, что может быть связано с затухающим интересом ко сладкому.\n",
    "\n",
    "Как и ожидалось, перед хэллоуином и новым годом (октябрь-декабрь) количество производимых сладостей существенно выше, чем в остальные месяцы, что вызвано повышенным спросом. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Выбор метрики"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Типичные метрики в задаче регрессии:\n",
    "- MAE\n",
    "- RMSE\n",
    "- MAPE\n",
    "- SMAPE\n",
    "\n",
    "Последние 2 метрики позволяют уйти от абсолютных величин к процентам, что является существенным преимуществом, когда целевая переменная меняется в большом диапазоне, как в решаемой задаче. SMAPE, вопреки своему названию, асимметрично штрафует за недопредсказания и перепредсказания, поэтому для оценки модели будем использовать MAPE.\n",
    "\n",
    "В обсуждениях XGBoost на github разработчики уточнили, что алгоритм плохо сходится при выборе MAPE в качестве метрики, поэтому также будем рассчитывать и MAE для дальнейшего сравнения моделей."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_absolute_percentage_error(y_true, y_pred): \n",
    "    y_true, y_pred = np.array(y_true), np.array(y_pred)\n",
    "    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Выбор модели"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Наиболее часто для предсказания временных рядов используются следующие методы:\n",
    "- Простейшие модели (наивная, предсказание средним, экспоненциальное сглаживание и т.п.)\n",
    "- Линейные модели\n",
    "- Эконометрические модели (ARIMA)\n",
    "- Facebook Prophet\n",
    "- XGBoost\n",
    "- LSTM RNNs\n",
    "\n",
    "Простейшие модели, как правило, используются как baseline и дают куда менее точные прогнозы, чем остальные методы.\n",
    "\n",
    "Линейные модели плохо справляются с нелинейными зависимостями и едва ли дадут хороший результат в решаемой задаче. Разве что, можно воспользоваться регрессией с lasso-регуляризагией для отбора признаков на вход другой модели.\n",
    "\n",
    "Чаще всего, в статьях на arxiv.org рассматриваемая задача решается с помощью LSTM-нейросеток, но мне пока что не хватает навыков в работе с нейросетями, поэтому данный метод также останется за рамками этой работы.\n",
    "\n",
    "А вот оставшиеся 3 модели можно проверить на имеющихся данных и выбрать ту, которая себя проявит лучше других."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Предобработка данных и создание новых признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ARIMA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сперва подготовим данные для обучения <b>ARIMA</b>-модели. В основе данного подхода лежит предположение о том, что исследуемый временной ряд стационарен. Если это не так, нужно произвести соответсвующие преобразования, чтобы он стал таковым."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Можно провести STL-декомпозицию ряда и проверить остатки на стационарность с помощью критерия Дики-Фуллера."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sm.tsa.seasonal_decompose(df['value']).plot()\n",
    "print(\"Критерий Дики-Фуллера: p=%f\" % sm.tsa.stattools.adfuller(df['value'])[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Критерий Дики-Фуллера отверг гипотезу о стационарности ряда. Следовательно, требуется его дальнейшая обработка."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Продиффериенциируем ряд с лагом 12 (по количеству месяцев)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['value_diff'] = df['value'] - df['value'].shift(12)\n",
    "sm.tsa.seasonal_decompose(df['value_diff'][12:]).plot()\n",
    "print(\"Критерий Дики-Фуллера: p=%f\" % sm.tsa.stattools.adfuller(df['value_diff'][12:])[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Остатки стационарны, тренд постоянен, можно подавать его на вход <b>ARIMA</b>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Отложим 2 года для дальнейшей оценки модели."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = df[:-24]\n",
    "df_test = df[-24:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Facebook Prophet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Главное преимущество <b>Facebook Prophet</b> - простота использования. Временной ряд не нуждается в предобработке и подаётся на вход в виде Pandas DataFrame'а с 2 колонками: 'ds' типа Date/Datetime и 'y' - значением прогнозируемой величины."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prophetdf = pd.DataFrame()\n",
    "prophetdf['ds'] = df.index\n",
    "prophetdf['y'] = df['value'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Также отложим 2 года для оценки модели."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prophetdf_train = prophetdf[:-24]\n",
    "prophetdf_test = prophetdf[-24:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XGBoost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для <b>XBGoost</b> нужно постараться нагенерировать побольше информативных категориальных и числовых признаков. Например, можно идти по временному ряду окном и считать различные статистики. Вместо того, чтобы делать это самому, а потом выбирать наиболее релевантные фичи из десятков а то и сотен, можно воспользоваться библиотекой <b>tsfresh</b>, которая делает это сама."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XGBdf, XGBy = make_forecasting_frame(df['value'], kind='kind', max_timeshift=36, rolling_direction=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XGBdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XGBX = extract_features(XGBdf, column_id=\"id\", column_sort=\"time\", column_value=\"value\", impute_function=impute,\n",
    "                     show_warnings=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(XGBX.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Отбросим константные признаки."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop constant features\n",
    "XGBX = XGBX.loc[:, XGBX.apply(pd.Series.nunique) != 1] \n",
    "print(XGBX.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Добавим предыдущее значение как признак."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XGBX[\"feature_last_value\"] = XGBy.shift(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XGBX = XGBX.iloc[1:, ]\n",
    "XGBy = XGBy.iloc[1: ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XGBX.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Мы получили 377 числовых признаков, характеризующих временной ряд, на которых уже можно обучить <b>XGBoost</b>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Снова отложим последние 2 года для оценки модели."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XGBX_train = XGBX[:-24]\n",
    "XGBy_train = XGBy[:-24]\n",
    "XGBX_test = XGBX[-24:]\n",
    "XGBy_test = XGBy[-24:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Кросс-валидация, настройка гиперпараметров модели, построение кривых валидации и обучения"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ARIMA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Очевидно, имеет смысл учесть в модели годовую сезонность - взять разность с лагом 12."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.subplot(211)\n",
    "sm.graphics.tsa.plot_acf(df_train['value_diff'][12:].values.squeeze(), lags=58, ax=ax)\n",
    "ax = plt.subplot(212)\n",
    "sm.graphics.tsa.plot_pacf(df_train['value_diff'][12:].values.squeeze(), lags=58, ax=ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Функцию подбора параметров модели можно позаимствовать из лекции. Границы интервалов поиска параметров ARIMA зададим исходя из коррелограмм, на которых явно прослеживается автокорреляция с лагом 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "results = []\n",
    "best_aic = float(\"inf\")\n",
    "\n",
    "ps = range(0, 3)\n",
    "d=1\n",
    "qs = range(0, 3)\n",
    "Ps = range(0, 3)\n",
    "D=1\n",
    "Qs = range(0, 3)\n",
    "parameters = product(ps, qs, Ps, Qs)\n",
    "parameters_list = list(parameters)\n",
    "\n",
    "for param in parameters_list:\n",
    "    try:\n",
    "        model=sm.tsa.statespace.SARIMAX(df_train['value'], order=(param[0], d, param[1]), \n",
    "                                        seasonal_order=(param[2], D, param[3], 12)).fit(disp=-1)\n",
    "    except ValueError:\n",
    "        print('wrong parameters:', param)\n",
    "        continue\n",
    "    aic = model.aic\n",
    "    if aic < best_aic:\n",
    "        best_model = model\n",
    "        best_aic = aic\n",
    "        best_param = param\n",
    "    results.append([param, model.aic])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_table = pd.DataFrame(results)\n",
    "result_table.columns = ['parameters', 'aic']\n",
    "print(result_table.sort_values(by = 'aic', ascending=True).head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Остатки лучшей модели:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(211)\n",
    "best_model.resid[13:].plot()\n",
    "plt.ylabel(u'Остатки')\n",
    "\n",
    "ax = plt.subplot(212)\n",
    "sm.graphics.tsa.plot_acf(best_model.resid[13:].values.squeeze(), lags=48, ax=ax)\n",
    "\n",
    "print(\"Критерий Стьюдента: p=%f\" % stats.ttest_1samp(best_model.resid[13:], 0)[1])\n",
    "print(\"Критерий Дики-Фуллера: p=%f\" % sm.tsa.stattools.adfuller(best_model.resid[13:])[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Исходя из критерия Стьюдента, остатки смещены. Однако они стационарны (подтверждается критерием Дики-Фуллера и визуально) и неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой). Неидеальный случай, но жить можно."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим, насколько хорошо модель описывает временной ряд."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train['model'] = best_model.fittedvalues\n",
    "df_train['value'].plot()\n",
    "df_train['model'][13:].plot(color='r')\n",
    "plt.ylabel('US candy production');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим, что модель очень хорошо подогналась под наши данные. Осталось проверить её на отложенной выборке."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Facebook Prophet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как уже упоминалось ранее, главное преимущество fbprophet - простота в использовании. Модель не требует настройки гиперпараметров."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fbmodel = fbprophet.Prophet(mcmc_samples=500)\n",
    "fbmodel.fit(prophetdf_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# XGBoost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обучим модель из не более чем 1000 итераций. Валидируемся разбиением тренировочной выборки на 5 фолдов."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtrain = xgb.DMatrix(XGBX_train, label=XGBy_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {\n",
    "        'objective': 'reg:linear',\n",
    "        'booster':'gblinear'\n",
    "    }\n",
    "rounds = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv = xgb.cv(params, dtrain, metrics = ('mae'), verbose_eval=False, nfold=5, \n",
    "            show_stdv=False, num_boost_round=rounds, seed=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bst = xgb.train(params, dtrain, num_boost_round=cv['test-mae-mean'].argmin())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим MAE на кросс-валидации."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv['test-mae-mean'].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction_train = bst.predict(dtrain)\n",
    "plt.plot(prediction_train)\n",
    "plt.plot(XGBy_train.values)\n",
    "plt.axis('tight')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Прогноз на отложенной выборке"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SARIMAX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train[\"arima_model\"] = best_model.fittedvalues\n",
    "forecast = best_model.predict(start = df_train.shape[0], end = df_train.shape[0]+24)\n",
    "forecast = df_train.arima_model.append(forecast).values\n",
    "forecast = forecast[12:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 7))\n",
    "plt.plot(forecast, color='r', label=\"model\")\n",
    "\n",
    "plt.title(\"SARIMA model\\n\")\n",
    "plt.plot(df['value'].values[12:], label=\"actual\")\n",
    "plt.legend()\n",
    "plt.axvspan(len(df['value'])-36, len(forecast), alpha=0.5, color='lightgrey')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видим, модель обманулась скачком значений прямо на конце тренировочного периода и на тесте выдала немного завышенный прогноз."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_absolute_percentage_error(df['value'].values[-24:], forecast[-24:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_absolute_error(df['value'].values[-24:], forecast[-24:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MAPE 10.00%, MAE 11.04 на отложенной выборке - довольно неплохой результат."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Facebook Prophet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "future = fbmodel.make_future_dataframe(periods=24, freq='M')\n",
    "fbforecast = fbmodel.predict(future)\n",
    "fbmodel.plot(fbforecast);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 7))\n",
    "plt.plot(fbforecast.yhat, color='r', label=\"model\")\n",
    "\n",
    "plt.title(\"fbprophet\\n\")\n",
    "plt.plot(prophetdf['y'].values[12:], label=\"actual\")\n",
    "plt.legend()\n",
    "plt.axvspan(len(df['value'])-36, len(fbforecast), alpha=0.5, color='lightgrey')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_absolute_percentage_error(prophetdf.y[-24:], fbforecast.yhat[-24:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_absolute_error(prophetdf.y[-24:], fbforecast.yhat[-24:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как итог, получаем MAPE в 11.55% и MAE в 12.46 на отложенной выборке. Как видно из графика, модель уловила основной тренд и годовую сезонность, но недостаточно резко реагирует на сильные колебания. Очень неплохой результат, учитывая, какими незначительными усилиями был получен прогноз, но другие модели позволяют достичь лучшего результата за счёт более тонкой настройки."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XGBoost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtest = xgb.DMatrix(XGBX_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction_test = bst.predict(dtest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction = np.append(prediction_train, prediction_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 7))\n",
    "plt.plot(prediction, color='r', label=\"model\")\n",
    "plt.title(\"XGBoost\\n\")\n",
    "plt.plot(df['value'].values[1:], label=\"actual\")\n",
    "plt.legend()\n",
    "plt.axvspan(len(df['value'])-36, len(fbforecast), alpha=0.5, color='lightgrey')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_absolute_error(XGBy_test, prediction_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "XGBoost наголову обогнал другие методы со вдвое меньшим значением MAE на отложенной выборке."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Выводы "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В ходе выполнения работы было получено 3 модели, способные делать прогноз уровня производства конфет в США на 2 года вперёд. Как оказалось, довольно неплохой прогноз можно получить почти без усилий при помощи библотки Facebook Prophet, а вот ARIMA требует дополнительной предобработки данных и настройки гиперпараметров. Ну а XGBoost вновь продемонстрировал, за что его так любят на Kaggle.\n",
    "\n",
    "В дальнейшем, прогнозы можно улучшить за счёт использования новой модели (например, LSTM RNN) или генерации других фич, используемых в обработке сигналов (например, Fast Fourier Transform)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
